Number and length of attractors in a critical Kauffman model with connectivity one 
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The Kauffman model describes a system of randomly connected nodes with dynamics based 
on Boolean update functions. Though it is a simple model it exhibits very complex behavior for 
"critical" parameter values at the boundary between a frozen and a disordered phase, and is therefore 
used for studies of real network problems. We prove here that the mean number and mean length of 
attractors in critical random Boolean networks with connectivity one both increase faster than any 
power law with network size. We derive these results by generating the networks through a growth 
process and by calculating lower bounds. 

PACS numbers: 89.75.Hc, 05.65.+b, 02.50.Cw 
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Boolean networks are often used as generic models for 
the dynamics of complex systems of interacting entities, 
such as social and economic networks, neural networks, 
and gene or protein interaction networks [Jj. The sim- 
plest and most widely studied of these models was intro- 
duced in 1969 by Kauffman as a model for gene reg- 
ulation. The system consists of N nodes, each of which 
receives input from K randomly chosen other nodes. The 
network is updated synchronously, the state of a node at 
time step t being a Boolean function of the states of the K 
input nodes at the previous time step, t—1. The Boolean 
updating functions are randomly assigned to every node 
in the network, and together with the connectivity pat- 
tern they define the realization of the network. For any 
initial condition, the network eventually settles on a pe- 
riodic attractor. Thus the number and the lengths of 
the attractors are important features of the networks. 
Of special interest are critical networks, which lie at the 
boundary between a frozen phase and a chaotic phase 
0, Q. In the frozen phase, a perturbation at one node 
propagates during one time step on an average to less 
than one node, and the attractor lengths remain finite in 
the limit N — > oo. In the chaotic phase, the difference be- 
tween two almost identical states increases exponentially 
fast, because a perturbation propagates on an average 
to more than one node during one time step Based 
on computer simulations, the mean attractor number of 
critical K — 2 Kauffman networks with a constant prob- 
ability distribution for the 16 possible updating functions 
was once believed to scale as y/~N Q. With increasing 
computer power, a faster increase was seen (linear in |6j, 
"faster than linear" in 0, stretched exponential in 
Finally, in a beautiful analytical study, Samuelsson and 
Troein |lfj| have proven that the number of attractors 
grows indeed faster than any power law with the network 
size N. Concerning the scaling behavior of the mean at- 
tractor length, there is not yet a conclusive result in the 
literature. While it appeared to increase as ^/N in earlier 
times |3, lUI , Bastolla and Parisi suggest that it might 
in fact increase faster than any power law [^, Q , and a 
recent review article treats this as an open question ||. 



Just as for the attractor number, computer simulations 
are hampered by undersampling, which makes it virtu- 
ally impossible to find attractors that occur only in few 
realizations or that have a small basin of attraction. 

In this letter, we want to prove that for K = 1 critical 
networks the mean number of attractors as well as their 
mean length grows faster than any power law with the 
network size. In such a network out of the four possible 
updating functions only the two non-constant ones occur. 
These networks are critical because a perturbation at one 
node propagates during one time step on an average to 
one node. This is the first analytical demonstration that 
in the same ensemble of networks both these quantities 
increase faster than any power law. If, as widely be- 
lieved, all critical Boolean networks behave in a similar 
way, these results should also hold for the critical K = 2 
networks. They can in fact directly be applied to the 
B\ class of the K — 2 critical networks, that is to the 
K = 2 networks with only those four Boolean functions 
that are canalizing and depend only on one of the two 
inputs. These networks are equivalent to critical K = 1 
networks since for Bi networks the effective number of in- 
puts per node is 1, as the other input does not influence 
the network dynamics. 

The topology of networks with K = 1 consists of loops 
and of trees rooted in them, and several exact results 
for these networks have been obtained by Flyvbjerg and 
Kjaer [1^, The dynamics on the loops determines the 
dynamics on the entire network, and the dynamics on 
the trees is slaved to the dynamics on the loops. Rele- 
vant nodes are those nodes that are not frozen and that 
control at least one other relevant element [|J. If all 4 
update functions are chosen with a nonzero probability, 
only short loops have a non-vanishing probability of not 
containing a constant function |l2|. Thus the number 
of relevant elements remains finite in the limit of infi- 
nite network size, and these networks are always in the 
frozen phase 0. Choosing only non-constant Boolean 
functions in K = 1 networks makes all nodes on all loops 
relevant, and all possible states on a loop are part of a 
cycle in state space. There are two kinds of non-constant 



2 



Boolean functions of one Boolean element: tautology (© 
coupling) and contradiction (G-coupling). The number 
and the lengths of the cycles on a loop depend only on 
the parity of the number of 0-functions and not on the 
details of the distribution of the Boolean functions. The 
cycle lengths of an "even" loop of length I are 1, I, and 
divisors of I. The maximum cycle length of an "odd" 
loop is 21 [12]. 

Let us first show that the mean attractor number in- 
creases faster than any power law with N. Let ni be the 
number of loops of length I, and m — Y]fLi ni I the num- 
ber of nodes in the loops, m is related to the attractors 



(1) 



where v\ is the number of attractors of length Ai . From 
here we see that finding an upper bound for the attractor 
length gives us a lower bound for the attractor number. 

The attractor length A is the least common multiple 
of cycle lengths (periods) of the loops, 
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For a fixed m, this product reaches its maximum if all 
U are equal, U — ZV«. In this case we have m = nil and 
A < 2 £™< . Maximizing this product as a function of the 
number ni of loops of length I 
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we obtain 

Am, 



<2]Jk < 2exp(— ) =2-2 
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A slight ly b etter upper bound of the form 2°- 5m was de- 
rived in |l2j, using a much more complicated calculation. 
From Eqs. © and QJ, we obtain a lower bound for the 
number Vi of attractors, 



(4) 



Averaging over the different network realizations gives 
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\^ v . > ± . 2 0A7m > ±2 0ATm . 
4^ ~ 2 ~ 2 
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An analytical expression for m can be derived from the 
exact results in [12J . One of them is the expectation value 
for the mean number of the loops of the length I in the 
large N limit 
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This result follows also from our Eqs. (|HJ and J7J below. 
Thus we find for the mean value of the number of nodes 
that are on loops 
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Approximating this sum with an integral one obtains 
m w y/^ N. Inserting this expression in Eq. J5|, we 
see that the mean number of attractors grows at least as 
fast as exp with the number of nodes. 

Next, we show that the mean attractor length diverges 
faster than any power law. For this purpose, we gener- 
ate the ensemble of all realizations of networks of size 
N + 1 via a growth process from the ensemble of net- 
works of size N. The following rule ensures that each 
network of the new ensemble is generated exactly once. 
The nodes are distinguishable and numbered in the se- 
quence in which they were added. To every network of 
the initial ensemble we insert a new node and add a new 
link. The new node has either itself as input or is linked 
to a node from the already existing network. Next, we 
have to assign to this new node all possible combinations 
of outgoing links. This is done such that all possible com- 
binations of the predecessor's outgoing links can become 
the outgoing links of the new node. If the new node is 
connected to itself any combination of outgoing links of 
node number 1 that are not on a loop can be moved to 
it. Different rearrangements of these links are weighed 
equally and every such network belongs to the new en- 
semble of the networks with N + 1 nodes. This procedure 
guarantees that the number of inputs per node of the al- 
ready existing network is not changed. If the node being 
the one linked to the inserted node was on the loop, there 
is a probability of \ that the new node is going to be on 
the loop (since \ is the probability that the predecessor's 
outgoing link that was the part of the loop is shifted to 
the new node). One can see that the already existing 
loops become bigger with time, and that new loops with 
only one node are created. We now consider the growth 
of the networks as a dynamical process, and we focus 
only on the loops. Since every node has the same dis- 
tribution of numbers and sizes of trees connected to it, 
all nodes in all loops become connected to a new node 
with the same probability in the ensemble. We define the 
time scale such that the rate of insertion of new nodes at 
a given position in a given loop is unity, 



dt = -^dN => t = hxVN. 
2N 



(7) 



Note that N now denotes the mean network size in the 
growing ensemble. By going from exact insertion num- 
bers to insertion rates, we have made a transition to a 
"grand canonical" ensemble. Within this ensemble, a 
loop of size I becomes a loop of size I + 1 with proba- 
bility Idt during a time interval dt. We then obtain the 
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following equations for the mean number n; of loops of 
size / 



df 



n\ = 1 — Hi 



ni = (I — — I ni \f I > 1 



With the initial condition A^(0) = 1, these equations have 
the solution 

HT(t) = 1 

Mt) = \-^ 2t +r~ 3t 

»H(t) = \ - | e -» + 2e- 3t -|e- 4t 



For the limiting case of large times, this solution can be 
approximated by 



1 I ~ 1 2/ 

WW = j - — e"» 
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which approaches the stationary solution nj = l^ 1 for 
t — > oo. Introducing the small parameter e « 1 as a 
measure of how far we are from the stationary solution, 
we find from 



lc 



that the critical value for the loop size for large eN is 



L = V2dV. 



(9) 



In the same manner we may write the master equation 
for the probability distribution P(n\, . . . , n{) of the loops 
smaller than I, 



— P(m, ...,ni) = ~y^i rij P(ni, . . .,m,rii + i, 

i—l 

I 

+ - l)(nj_i + 1) P(ni, . . .,n,_i + l,n i+ i, . . . ,nj) 

i=l 

The stationary solution for this expression valid for the 
loops smaller than l c is 

l / 1 \ ™< i 

i— 1 x 7 

This solution is time independent and we can conclude 
that the distribution of the loops smaller than l c is not 
changing with time, i.e. with the growth of the system 



size. Furthermore the probabilities for having rij loops 
of size i are independent from each other and Poisson 
distributed with a mean i . 

Equipped with these results, we can now evaluate the 
lower bound for the mean attractor length A. Suppose 
that the system is enlarged so that its number of nodes 
is N' = aN, with a > 1 and N nodes of the previous 
system. The length of the attractor is the least com- 
mon multiple of the loop periods, i.e. the cycle lengths 
of the loops. Since Ai 
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we obtain a 



lower bound by evaluating only the change of the least 
common multiple of the periods of loops smaller than 

with increas- 
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l c , that is the change of A; i <i c _ ^ , 2e}f , 
ing system size. Our above considerations show that the 
distribution of loops of size smaller than l c (N) does not 
change when going to an ensemble of systems of size N'. 
However, these systems contain additional loops in the 
interval [l c (N), l c (N')]. If the period of such an addi- 
tional loop is a prime number larger than l c (N), the least 
common multiple of all loop periods is multiplied by this 
period. A loop with a prime number of nodes has only 
two possible periods: 1 or I if the loop is even, and 2 
or 21 if the loop is odd. If the additional loop is not on 
the cycle of length 1 or 2, the least common multiple of 
the periods of the loops smaller than l c (N') is at least as 
large as the product of the new loop size, and the least 
common multiple of the periods of the loops smaller than 
l c (N). The number of primes not exceeding the value of 
some positive number x is asymptotically expressed as 
n x = x/lnx (see, e.g. ^^)- The probability that a ran- 
domly chosen number in the interval [l c (N),l c (N')j is a 
prime number is 
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This is identical to the probability that the new loop 
size is a prime number. Taking all these considerations 
together, we have 
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The probability Pi 00 p for having a loop in the interval 
[l c (N),l e (N')] is obtained using (TOjl. The probability of 
having no loop of size I, m = 0, is e' 1 ' 1 . Thus the upper 
bound for the probability of having no loops with size 
from the interval [l c (N),l c (N')] is 
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The probability P ri oti,2 that the new loop is not on an 
attractor of length 1 or 2 is obtained as follows: The 
number of its cycles is (2 Z — 2)/Z + 2 in the case of the even 
loops and (2 l — 2)/2l+ 1 for the odd loops. Among these 
cycles two are of length 1 for the first type of loop and 
one is of length 2 for the second type. The probability 
that the loop of size I is not on a cycle of length 1 or 2 is 
1 — 21 /2 l for large values of /. The loop we arc observing 
is of size y/2eN and for the probability P no ti,2 we obtain 



2V27N 
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(13) 



For a given e this probability is non- vanishing, 
i.e. P no ti,2 > V > 0, if N > 2/e. Since we are considering 
the limit of large N, this condition is satisfied. Applying 
this result to the lower bound for the attractor size we 
finally have 
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Setting N = a^No and defining a constant C 
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Eq. I|14|) can be transformed into 
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and finally with \i = (\n(N / N ) / \n(a)) into 
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This increases faster than any power law with N, but 
slower than a stretched exponential. 

It remains to see in how far these results apply also to 
critical K — 2 Kauffman networks. Bastolla and Parisi 
■9f have pointed out that the set of relevant nodes of these 
networks consists of modules, which together determine 
the attractor numbers and lengths. This situation is very 
similar to the critical K — 1 networks treated in this pa- 
per, where the modules built of relevant nodes are loops, 
the properties of which determine the attractors. The 
main difference in K = 2 networks is that the modules 



are more complicated. A thorough evaluation of their 
properties has not yet been done. 

Nevertheless, the evidence cited so far suggests that 
for all critical Kauffman networks the mean number and 
length of attractors diverges faster than any power law. 
This means that the attractors are too long and too many 
to represent cellular differentiation, to which the model 
was originally applied. The vast number of attractors 
in these models appears to be a consequence of the syn- 
chronous updating scheme. Recent studies of modified 
models that allow for randomness in the updating rules, 
indicate that a deviation from synchronous update re- 
duces the number of attractors considerably [14(, which 
now becomes a power law |l5l llfij . However, in order to 
model biological networks realistically, further modifica- 
tions are needed, and the present work is only a small 
step on the long way towards understanding regulatory 
networks. 

We thank Viktor Kaufman for useful discussions. 
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